home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Analog / LSolve.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  10.7 KB  |  335 lines
  1. (*  :Title:    Solution of Differential Equations via Laplace ransforms. *)
  2.  
  3. (*  :Authors:    Wally McClure, Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    This package will solve linear constant coefficient
  7.         differential equations by means of the Laplace transform.
  8.  *)
  9.  
  10. (*  :Context:    SignalProcessing`Analog`LSolve`  *)
  11.  
  12. (*  :PackageVersion:  2.7    *)
  13.  
  14. (*
  15.     :Copyright:    Copyright 1990-1991 by Brian L. Evans
  16.         Georgia Tech Research Corporation
  17.  
  18.     Permission to use, copy, modify, and distribute this software
  19.     and its documentation for any purpose and without fee is
  20.     hereby granted, provided that the above copyright notice
  21.     appear in all copies and that both that copyright notice and
  22.     this permission notice appear in supporting documentation,
  23.     and that the name of the Georgia Tech Research Corporation,
  24.     Georgia Tech, or Georgia Institute of Technology not be used
  25.     in advertising or publicity pertaining to distribution of the
  26.     software without specific, written prior permission.  Georgia
  27.     Tech makes no representations about the suitability of this
  28.     software for any purpose.  It is provided "as is" without
  29.     express or implied warranty.
  30.  *)
  31.  
  32. (*
  33.     :History:    Date started -- September 28, 1990 (adapted from ZSolve.m)
  34.         Handle arbitrary initial conditions -- May, 1991
  35.  *)
  36.  
  37. (*  :Keywords:    linear constant-coefficient differential equations    *)
  38.  
  39. (*  :Source:    {Operational Mathematics} by Ruel Churchill  *)
  40.  
  41. (*  :Warning:    *)
  42.  
  43. (*  :Mathematica Version:  1.2 or 2.0  *)
  44.  
  45. (*  :Limitation:  *)
  46.  
  47. (*  :Discussion:  *)
  48.  
  49. (*  :Functions:   PartialL  LSolve  *)
  50.  
  51.  
  52.  
  53. (*  B E G I N     P A C K A G E  *)
  54.  
  55. BeginPackage[ "SignalProcessing`Analog`LSolve`",
  56.           "SignalProcessing`Analog`LaPlace`",
  57.           "SignalProcessing`Analog`InvLaPlace`",
  58.           "SignalProcessing`Analog`LSupport`",
  59.           "SignalProcessing`Support`TransSupport`",
  60.           "SignalProcessing`Support`ROC`",
  61.           "SignalProcessing`Support`SigProc`",
  62.           "SignalProcessing`Support`SupCode`" ]
  63.  
  64.  
  65. If [ TrueQ[ $VersionNumber >= 2.0 ],
  66.      Off[ General::spell ];
  67.      Off[ General::spell1 ] ];
  68.  
  69.  
  70. (*  U S A G E     I N F O R M A T I O N  *)
  71.  
  72. LSolve::usage =
  73.     "LSolve[ diffequ == drivingfun, y[t] ] solves the differential \
  74.     equation diffequ = drivingfun, where diffequ is a linear constant \
  75.     coefficient differential equation and drivingfun is the driving \
  76.     function (a function of t). \
  77.     Thus, diffequ has the form a0 y[t] + a1 y'[t] + .... \
  78.     One can specify initial values; e.g., \
  79.     LSolve[ y''[t] + 3/2 y'[t] + 1/2 y[t] == Exp[a t], \
  80.     y[t], y[0] -> 4, y'[0] -> 10 ]. \
  81.     A differential equation of N terms needs N-1 initial conditions. \
  82.     All unspecified conditions are considered to be zero. \
  83.     LSolve can justify its answers."
  84.  
  85. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  86.  
  87.  
  88. Begin[ "`Private`" ]
  89.  
  90.  
  91. (*  Messages  *)
  92. LSolve::nolaplace =
  93.     "The forcing function `` does not have a valid Laplace ransform."
  94. LSolve::toomany =
  95.     "Initial/boundary conditions occur at more than two locations in ``."
  96. PartialL::nolaplace =
  97.     "The differential equation `` does not have a partial Laplace transform."
  98.  
  99.  
  100. (*  completetransform  *)
  101. completetransform[ a_. y_ + b_., y_, x_ ] := ( x - b ) / a /; FreeQ[b, y]
  102.  
  103.  
  104. (* ic2index --  extract location of initial/bondary condition *)
  105. SetAttributes[ic2index, Listable]
  106. ic2index[ (y_[index_] -> value_), y_ ] := index
  107. ic2index[ (Derivative[n_][y_][index_] -> value_), y_ ] :=  index
  108.  
  109.  
  110. (* introduction *)
  111. introduction[dialogue_, diffequ_, drivingfun_, bvflag_, initlist_, y_[t_]] :=
  112.     Block [ {},
  113.         If [ dialogue ||  SameQ[dialogue, All],
  114.              Print["Solving for ",y[t]," in the differential equation"];
  115.              Print["  ", diffequ, " = ", drivingfun];
  116.              If [ bvflag,
  117.                   Print["  subject to the boundary conditions"],
  118.                   Print["  subject to the initial conditions"] ];
  119.              If [ EmptyQ[initlist],
  120.               Print["  being zero."],
  121.                   Print["  ", initlist ] ] ] ]
  122.  
  123.  
  124. (* preprocess *)
  125. preprocess[ left_ == right_, y_[t_] ] :=
  126.     Block [    {diffequ, drivingfun},
  127.  
  128.         diffequ = Select[ left - right, validterm[y[t]] ];
  129.         drivingfun = -(left - right - diffequ);
  130.  
  131.         diffequ == drivingfun ]
  132.  
  133.  
  134. (* printdialogue *)
  135. printdialogue[dialogue_, y_[t_], Ypartial_, X_, s_, bvflag_] :=
  136.     Block [    {formula, ll, n, rules, term1, term2, term3, terml, ul, YYstr},
  137.  
  138.         (* formula for Laplace transform of nth derivative of y(t) *)
  139.         YYstr = ToString[StringForm["L { `` }", y[t]]];
  140.         term1 = s^n YYstr;
  141.         If [ bvflag,
  142.              rules = { ul -> "ul", ll -> "ll" };
  143.              term2 = s^(n-1) ( y[ul] Exp[-s ul] - y[ll] Exp[-s ll] );
  144.              term3 = s^(n-2) ( y'[ul] Exp[-s ul] - y'[ll] Exp[-s ll] );
  145.              terml = ( Derivative[n-1][y][ul] Exp[-s ul] -
  146.                    Derivative[n-1][y][ll] Exp[-s ll] ),
  147.  
  148.              rules = { ll -> "t0" };
  149.              term2 = - s^(n-1) y[ll] Exp[-s ll];
  150.              term3 = - s^(n-2) y'[ll] Exp[-s ll];
  151.              terml = - Derivative[n-1][y][ll] Exp[-s ll] ];
  152.         rules = rules ~Join~ { s -> "s", n -> "n" };
  153.  
  154.         term1 = term1 /. rules;
  155.         term2 = term2 /. rules;
  156.         term3 = term3 /. rules;
  157.         terml = terml /. rules;         (* yeh, this was a lot of work *)
  158.  
  159.         (* dialogue for user *)
  160.         If [ dialogue,
  161.              Print["After taking the Laplace transform of both sides"];
  162.              Print["and solving for the transform of ", y[t], ":" ] ];
  163.         If [ SameQ[dialogue, All],
  164.              Print[ " " ];
  165.              Print[ "The Laplace transform of the left side is:" ];
  166.              Print[ "  ",
  167.                     ( Ypartial /. YY -> YYstr ) /. s -> "s" ];
  168.              Print[ "( In the general case, the unilateral Laplace" ];
  169.              Print[ "  transform of the nth derivative of ",
  170.                 y[t], " is:" ];
  171.              Print[ "  L {", Derivative["n"][y][t], " } = ", term1,
  172.                 " + ", term2, " + ", term3, " + ... + ", terml ];
  173.              If [ bvflag,
  174.                   Print["  where t=ll is the lower limit and t=ul is"];
  175.                 Print["  the upper limit of the boundary conditions. )"],
  176.               Print["  where t=t0 is the initial condition. )"] ];
  177.              Print[ " " ];
  178.              Print[ "The Laplace transform of the right side is:" ];
  179.              Print[ "  ", X /. s -> "s" ];
  180.              Print[ "Solving for the unknown transform yields" ] ] ]
  181.  
  182.  
  183. (* validic --  valid initial/bondary condition? *)
  184. validic [y_[t_]] [ y_[index_] -> value_ ] := True
  185. validic [y_[t_]] [ Derivative[n_][y_][index_] -> value_ ] := True
  186.  
  187.  
  188. (* validterm --  valid term in left-hand side of a differential equation? *)
  189. validterm [y_[t_]] [ a_. y_[t_ + t0_.] ] := True
  190. validterm [y_[t_]] [ a_. Derivative[n_][y_][t_] ] := True
  191.  
  192.  
  193.  
  194.  
  195. (*  PartialL  *)
  196. PartialL[ x_ + v_, rest__ ] :=
  197.     PartialL[ x, rest ] + PartialL[ v, rest ]
  198.  
  199. PartialL[ a_ x_, y_[t_], s_, YY_, rest___ ] :=
  200.     a PartialL[ x, y[t], s, YY, rest ] /; FreeQ[a, t]
  201.  
  202. PartialL[ Derivative[1][y_][t_ + a_.], y_[t_], s_, YY_, ll_, Infinity ] :=
  203.     Exp[a s] ( s YY - y[ll] Exp[-s ll] ) /;
  204.     FreeQ[a, t]
  205.  
  206. PartialL[ Derivative[1][y_][t_ + a_.], y_[t_], s_, YY_, ll_, ul_ ] :=
  207.     Exp[a s] ( s YY + y[ul] Exp[-s ul] - y[ll] Exp[-s ll] ) /;
  208.     FreeQ[a, t] && ! SameQ[ul, Infinity]
  209.  
  210. PartialL[ Derivative[n_][y_][t_ + a_.], y_[t_], s_, YY_, ll_, Infinity ] :=
  211.     Block [    {i, otherterms, term, term1, term2},
  212.         term1 = s^n YY;
  213.         term2 = - s^(n-1) y[ll] Exp[-s ll];
  214.         term = - Derivative[i][y][ll] Exp[-s ll];
  215.         otherterms = Apply[Plus, Table[s^(n-i-1) term, {i, 1, n-1}]];
  216.         Exp[a s] ( term1 + term2 + otherterms ) ] /;
  217.     FreeQ[a, t] && n > 1
  218.  
  219. PartialL[ Derivative[n_][y_][t_ + a_.], y_[t_], s_, YY_, ll_, ul_ ] :=
  220.     Block [    {i, otherterms, term, term1, term2},
  221.         term1 = s^n YY;
  222.         term2 = s^(n-1) ( y[ul] Exp[-s ul] - y[ll] Exp[-s ll] );
  223.         term = ( Derivative[i][y][ul] Exp[-s ul] -
  224.              Derivative[i][y][ll] Exp[-s ll] );
  225.         otherterms = Apply[Plus, Table[s^(n-i-1) term, {i, 1, n-1}]];
  226.         Exp[a s] ( term1 + term2 + otherterms ) ] /;
  227.     FreeQ[a, t] && n > 1 && ! SameQ[ul, Infinity]
  228.  
  229. PartialL[ y_[t_ + a_.], y_[t_], s_, YY_, rest___ ] :=
  230.     Exp[a s] YY /; FreeQ[a, t]
  231.  
  232. PartialL[ a_, y_[t_], s_, YY_, rest___ ] :=
  233.     Message[ PartialL::nolaplace, a ] /; FreeQ[a, t]
  234.  
  235.  
  236.  
  237.  
  238. (*  LSolve  *)
  239. LSolve/: Options[LSolve] := { Dialogue -> True }
  240.  
  241. LSolve[ y_[t_ + t0_.] == drivingfun_, y_[t_] ] :=
  242.     ( drivingfun /. { t -> t - t0 } ) /; FreeQ[t0, t]
  243.  
  244. LSolve[ left_ == right_, y_[t_], init___ ] :=
  245.     Block [    {auginitlist, bvflag, dialogue, diffequ, drivingfun,
  246.          equationlist, indexlist, initfun, lowerlimit, options,
  247.          rminus, rplus, s, sobj, upperlimit, X, Y, Ycollected,
  248.          Ypartial, YY, YYstr},
  249.  
  250.         (* Parse options, determine the level of Dialogue requested *)
  251.         options = ToList[init] ~Join~ Options[LSolve];
  252.         dialogue = Replace[Dialogue, options];
  253.         initlist = Select[ options, validic[y[t]] ];
  254.         auginitlist = initlist ~Join~
  255.                   { y[t0_] :> 0, Derivative[n_][y][t0_] :> 0 };
  256.  
  257.         (* Collect terms  *)
  258.         equationlist = preprocess[left == right, y[t]];
  259.         diffequ = First[equationlist];
  260.         drivingfun = Second[equationlist];
  261.  
  262.         (* Find the location in t of initial/boundary conditions *)
  263.         If [ EmptyQ[initlist],
  264.              bvflag = False;
  265.              lowerlimit = 0;
  266.              upperlimit = Infinity,
  267.  
  268.              indexlist = Union[ ic2index[initlist, y] ];
  269.              If [ Length[indexlist] > 2,
  270.                   Message[ LSolve::toomany, t ];
  271.                   indexlist = Take[indexlist, 2] ];
  272.              indexlist = Sort[indexlist];
  273.              lowerlimit = First[indexlist];
  274.              bvflag = ( Length[indexlist] == 2 );
  275.              upperlimit = If [ bvflag, Second[indexlist], Infinity ] ];
  276.  
  277.         (* Tell the user about the differential equation *)
  278.         introduction[ dialogue, diffequ, drivingfun,
  279.                   bvflag, initlist, y[t] ];
  280.  
  281.         (* Find complete Laplace transform of driving function *)
  282.         sobj = LaPlace[ drivingfun, t, s, Dialogue -> False ];
  283.         If [ ! SameQ[Head[sobj], LTransData],
  284.              Message[ LSolve::nolaplace, drivingfun ];
  285.              Return[] ];
  286.         rminus = GetRMinus[sobj];              (* R- of X(s) *)
  287.         rplus = GetRPlus[sobj];                  (* R+ of X(s) *)
  288.         X = TheFunction[sobj];        (* X(s) from Laplace object *)
  289.  
  290.         (* Find partial Laplace transform of differential equation *)
  291.         Ypartial = PartialL[ diffequ, y[t], s, YY,
  292.                      lowerlimit, upperlimit ];
  293.         Ycollected = Collect[Ypartial, YY];
  294.  
  295.         (* Dialogue taking Laplace transform of both sides to user *)
  296.         printdialogue[dialogue, y[t], Ycollected, X, s, bvflag];
  297.  
  298.         (* Inverse transform complete Laplace transform of y, Y(s) *)
  299.         Y = completetransform[Ycollected, YY, X];
  300.  
  301.         If [ SameQ[dialogue, All],
  302.              Print[ "  ", Y /. s -> "s" ];
  303.              Print[ "which becomes" ] ];
  304.  
  305.         (*  Replace initial conditions in Y(s) with values *)
  306.         Y = Y /. auginitlist;
  307.  
  308.         If [ dialogue || SameQ[dialogue, All],
  309.              Print[ "  ", Y /. s -> "s" ];
  310.              Print[ "Inverse transforming this gives ", y[t], ":" ] ];
  311.  
  312.         InvLaPlace[{Y, rminus, rplus}, s, t, Dialogue -> False] ]
  313.  
  314.  
  315. (*  E N D     P A C K A G E  *)
  316.  
  317. End[]
  318. EndPackage[]
  319.  
  320. If [ TrueQ[ $VersionNumber >= 2.0 ],
  321.      On[ General::spell ];
  322.      On[ General::spell1 ] ];
  323.  
  324.  
  325. (*  H E L P     I N F O R M A T I O N  *)
  326.  
  327. Combine[ SPfunctions, { LSolve } ]
  328. Protect[ LSolve ]
  329.  
  330.  
  331. (*  E N D I N G     M E S S A G E  *)
  332.  
  333. Print[ "LSolve, a differential equation solver, is loaded." ]
  334. Null
  335.